En este ejercicio se va a estimar un modelo VAR en donde se relaciona la inflación y la depreciación para Colombia. Los datos son tomados de Posada y Londoño-Sierra (2022) y se cuenta con la siguiente información mensual de enero de 2003 a mayo de 2022 \((n=224)\):

Las variables endógenas son la inflación para Colombia y la depreciación, mientras como variable exógena se tiene la inflación de Estados Unidos.

Análisis gráfico

A continuación cargamos los datos y realizamos los gráficos de las series

library(fpp3); library(tidyverse); library(vars); library(urca); library(forecast)
library(openxlsx); library(plotly)
data <- read.xlsx("https://ebadilloe.github.io/EconometriaII/IntroVar/Inflacion_depreciacion_enero2003_junio2022.xlsx") |> 
  mutate(date = as.yearmon(date))

ggplot(data) +
  geom_line(aes(x = date, y = infcol, colour = "Colombia"), size = 0.8) +
  geom_line(aes(x = date, y = Infeeuu, colour = "Estados Unidos"), size = 0.8) +
  geom_line(aes(x = date, y = Dep/4, colour = "Depreciación"), size = 0.8) +
  scale_colour_manual(name = "", values = c("Colombia" = "darkblue", "Estados Unidos" = "red", "Depreciación" = "green"), limits = c("Colombia", "Estados Unidos", "Depreciación")) +
  labs(title = "Gráfico 1. Inflación en Colombia y Estados Unidos, y depreciación", x = "Mes", y ="") + 
  theme(legend.position = c(0.17, 0.93),
        legend.background = element_rect(fill="transparent", colour ="transparent"),
        legend.key = element_rect(fill = "transparent", colour = "transparent")) +
  scale_y_continuous(name = "Inflación (%)",
                     sec.axis = sec_axis(~.*4, name="Depreciación (%)"))

Es posible generar un gráfico interactivo para ver los valores directamente en el gráfico

g1 <- ggplot(data) +
  geom_line(aes(x = date, y = infcol, colour = "Colombia"), size = 0.8) +
  geom_line(aes(x = date, y = Infeeuu, colour = "Estados Unidos"), size = 0.8) +
  scale_colour_manual(name = "", values = c("Colombia" = "darkblue", "Estados Unidos" = "red"), limits = c("Colombia", "Estados Unidos")) +
  labs(title = "Gráfico 1. Inflación en Colombia y Estados Unidos, y depreciación", x = "Mes", y ="") + 
  theme(legend.position = c(0.17, 0.93),
        legend.background = element_rect(fill="transparent", colour ="transparent"),
        legend.key = element_rect(fill = "transparent", colour = "transparent")) +
  scale_y_continuous(name = "Inflación (%)")

ay <- list(
  tickfont = list(size=12),
  titlefont=list(size=15),
  overlaying = "y",
  nticks = 5,
  zeroline=F,
  side = "right",
  title = "Depreciación (%)"
)

ggplotly(g1) |> 
  add_lines(x=~date, y=~Dep, yaxis="y2", line = list(color = "#00FF00", width = 3), 
            data=data, showlegend=T, name = "Depreciación") |> 
  layout(yaxis2 = ay)

Pruebas de estacionariedad

Ahora pasamos a determinar si las series son o no estacionarias y para ello calculamos los correlogramas y las pruebas de raices unitarias

ggAcf(data$infcol, main="Función AC inflación Colombia", 
      ylab="AC")

ggPacf(data$infcol, main="Función PAC inflación Colombia", 
       ylab="PAC")

ggAcf(data$Infeeuu, 
      main="Función AC inflación Estados Unidos", 
      ylab="AC")

ggPacf(data$Infeeuu, 
       main="Función PAC inflación Estados Unidos", 
       ylab="PAC")

ggAcf(data$Dep, main="Función AC depreciación", 
      ylab="AC")

ggPacf(data$Dep, main="Función PAC depreciación", 
       ylab="PAC")

DF_infcol <- ur.df(data$infcol, type="trend", 
                   selectlags = "AIC")
summary(DF_infcol)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01529 -0.19441 -0.02086  0.15108  1.10567 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0420082  0.0740164   0.568    0.571    
## z.lag.1     -0.0162162  0.0113646  -1.427    0.155    
## tt           0.0003089  0.0003133   0.986    0.325    
## z.diff.lag   0.5362034  0.0578482   9.269   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2971 on 228 degrees of freedom
## Multiple R-squared:  0.2955, Adjusted R-squared:  0.2863 
## F-statistic: 31.88 on 3 and 228 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.4269 1.5218 2.2241 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
DF_infeeuu <- ur.df(data$Infeeuu, type="trend", 
                    selectlags = "AIC")
summary(DF_infeeuu)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.96687 -0.23152  0.01454  0.22984  1.47638 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0491637  0.0684954   0.718   0.4736    
## z.lag.1     -0.0412450  0.0172245  -2.395   0.0174 *  
## tt           0.0005097  0.0004149   1.229   0.2205    
## z.diff.lag   0.4869905  0.0602511   8.083  3.7e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4209 on 228 degrees of freedom
## Multiple R-squared:  0.2341, Adjusted R-squared:  0.2241 
## F-statistic: 23.23 on 3 and 228 DF,  p-value: 3.688e-13
## 
## 
## Value of test-statistic is: -2.3945 2.505 3.6245 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
DF_dep <- ur.df(data$Dep, type="trend", selectlags = "AIC")
summary(DF_dep)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.7472  -2.8904   0.0466   2.2666  19.2602 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.876071   0.702580  -1.247 0.213701    
## z.lag.1     -0.097497   0.025608  -3.807 0.000181 ***
## tt           0.009415   0.005380   1.750 0.081432 .  
## z.diff.lag   0.113791   0.064996   1.751 0.081338 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.228 on 228 degrees of freedom
## Multiple R-squared:  0.06642,    Adjusted R-squared:  0.05414 
## F-statistic: 5.407 on 3 and 228 DF,  p-value: 0.001304
## 
## 
## Value of test-statistic is: -3.8072 4.9807 7.4484 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

La conclusión es que las series no son estacionarias y es necesario diferenciarlas

data <- data |> 
  mutate(dinfcol = c(NA,diff(infcol,1)),
         dinfeeuu = c(NA,diff(Infeeuu,1)),
         ddep = c(NA,diff(Dep,1)))
plot.ts(data[,5:7])

Realizamos pruebas de raices unitarias a las series diferenciadas para chequear que ya sean estacionarias

DF_dinfcol <- ur.df(na.omit(data$dinfcol), type="trend", 
                   selectlags = "AIC")
summary(DF_dinfcol)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09106 -0.19032 -0.00871  0.16202  1.07029 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0549276  0.0399652  -1.374   0.1707    
## z.lag.1     -0.4781410  0.0649994  -7.356 3.43e-12 ***
## tt           0.0005113  0.0002987   1.712   0.0883 .  
## z.diff.lag  -0.0022002  0.0667830  -0.033   0.9737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2973 on 227 degrees of freedom
## Multiple R-squared:  0.2395, Adjusted R-squared:  0.2294 
## F-statistic: 23.83 on 3 and 227 DF,  p-value: 1.92e-13
## 
## 
## Value of test-statistic is: -7.3561 18.0765 27.1132 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
DF_dinfeeuu <- ur.df(na.omit(data$dinfeeuu), type="trend", 
                    selectlags = "AIC")
summary(DF_dinfeeuu)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01980 -0.21012  0.01721  0.20129  1.51427 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0477356  0.0551258  -0.866 0.387438    
## z.lag.1     -0.6819312  0.0674218 -10.114  < 2e-16 ***
## tt           0.0005704  0.0004112   1.387 0.166755    
## z.diff.lag   0.2468993  0.0642710   3.842 0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4135 on 227 degrees of freedom
## Multiple R-squared:  0.3167, Adjusted R-squared:  0.3076 
## F-statistic: 35.07 on 3 and 227 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.1144 34.11 51.159 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47
DF_ddep <- ur.df(na.omit(data$ddep), type="trend", selectlags = "AIC")
summary(DF_ddep)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.2610  -2.9848   0.3841   2.5277  17.9392 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.514409   0.718268  -0.716    0.475    
## z.lag.1     -0.963697   0.090482 -10.651   <2e-16 ***
## tt           0.003649   0.005333   0.684    0.494    
## z.diff.lag   0.037577   0.066754   0.563    0.574    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.397 on 227 degrees of freedom
## Multiple R-squared:  0.4654, Adjusted R-squared:  0.4583 
## F-statistic: 65.86 on 3 and 227 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.6507 37.8171 56.7256 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

Los tests rechazan la \(H_0\): existe raíz unitaria, lo que implica que las series diferenciadas son estacionarias, es decir, son \(I(1)\). Por lo tanto, en el modelo VAR las series se incorporan en diferencias

Estimación del modelo VAR

Ahora se determina el rezago óptimo del VAR

data.var <- cbind(na.omit(data$dinfcol), na.omit(data$ddep))
colnames(data.var) <- c("dinfcol", "ddep")
exo <- as.matrix(data.frame(na.omit(data$dinfeeuu)))
colnames(exo) <- "dinfeeuu"

select.var <- VARselect(data.var, 
                        lag.max = 20, 
                        type = "const",
                        exogen = exo)
select.var$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     13      1      1     13

El número de rezagos óptimo puede ser 1 o 13 así que se estiman los dos modelos.

var1.est <- VAR(data.var, p = 1, type = "const", exogen = exo)
summary(var1.est)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dinfcol, ddep 
## Deterministic variables: const 
## Sample size: 232 
## Log Likelihood: -748.511 
## Roots of the characteristic polynomial:
## 0.5123 0.06524
## Call:
## VAR(y = data.var, p = 1, type = "const", exogen = exo)
## 
## 
## Estimation results for equation dinfcol: 
## ======================================== 
## dinfcol = dinfcol.l1 + ddep.l1 + const + dinfeeuu 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## dinfcol.l1 0.536828   0.052966  10.135  < 2e-16 ***
## ddep.l1    0.016882   0.003569   4.731 3.93e-06 ***
## const      0.003424   0.018489   0.185    0.853    
## dinfeeuu   0.167732   0.040164   4.176 4.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.2811 on 228 degrees of freedom
## Multiple R-Squared: 0.3692,  Adjusted R-squared: 0.3609 
## F-statistic: 44.48 on 3 and 228 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation ddep: 
## ===================================== 
## ddep = dinfcol.l1 + ddep.l1 + const + dinfeeuu 
## 
##            Estimate Std. Error t value Pr(>|t|)  
## dinfcol.l1 -0.68427    1.00797  -0.679   0.4979  
## ddep.l1     0.04074    0.06791   0.600   0.5492  
## const      -0.03329    0.35185  -0.095   0.9247  
## dinfeeuu   -1.44274    0.76434  -1.888   0.0604 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.35 on 228 degrees of freedom
## Multiple R-Squared: 0.02234, Adjusted R-squared: 0.009473 
## F-statistic: 1.736 on 3 and 228 DF,  p-value: 0.1603 
## 
## 
## 
## Covariance matrix of residuals:
##         dinfcol    ddep
## dinfcol 0.07904  0.1027
## ddep    0.10267 28.6235
## 
## Correlation matrix of residuals:
##         dinfcol    ddep
## dinfcol 1.00000 0.06826
## ddep    0.06826 1.00000
var13.est <- VAR(data.var, p = 13, type = "const", exogen = exo)
summary(var13.est)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dinfcol, ddep 
## Deterministic variables: const 
## Sample size: 220 
## Log Likelihood: -631.007 
## Roots of the characteristic polynomial:
## 0.9714 0.9714 0.9632 0.9632 0.9626 0.9626 0.9367 0.9367 0.9355 0.9355 0.9274 0.9274 0.9252 0.9252 0.9219 0.9219 0.903 0.903 0.8968 0.8968 0.8844 0.8844 0.8805 0.8805 0.3378 0.278
## Call:
## VAR(y = data.var, p = 13, type = "const", exogen = exo)
## 
## 
## Estimation results for equation dinfcol: 
## ======================================== 
## dinfcol = dinfcol.l1 + ddep.l1 + dinfcol.l2 + ddep.l2 + dinfcol.l3 + ddep.l3 + dinfcol.l4 + ddep.l4 + dinfcol.l5 + ddep.l5 + dinfcol.l6 + ddep.l6 + dinfcol.l7 + ddep.l7 + dinfcol.l8 + ddep.l8 + dinfcol.l9 + ddep.l9 + dinfcol.l10 + ddep.l10 + dinfcol.l11 + ddep.l11 + dinfcol.l12 + ddep.l12 + dinfcol.l13 + ddep.l13 + const + dinfeeuu 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## dinfcol.l1   0.6373281  0.0685471   9.298  < 2e-16 ***
## ddep.l1      0.0125264  0.0038727   3.235 0.001435 ** 
## dinfcol.l2  -0.1285394  0.0764475  -1.681 0.094309 .  
## ddep.l2      0.0021218  0.0034600   0.613 0.540449    
## dinfcol.l3   0.1826525  0.0749046   2.438 0.015659 *  
## ddep.l3      0.0021072  0.0034400   0.613 0.540886    
## dinfcol.l4  -0.0915913  0.0748387  -1.224 0.222509    
## ddep.l4     -0.0025055  0.0034213  -0.732 0.464860    
## dinfcol.l5   0.0619781  0.0752342   0.824 0.411074    
## ddep.l5      0.0060039  0.0034459   1.742 0.083050 .  
## dinfcol.l6   0.0901208  0.0748853   1.203 0.230283    
## ddep.l6     -0.0048866  0.0033989  -1.438 0.152140    
## dinfcol.l7  -0.0723491  0.0747549  -0.968 0.334353    
## ddep.l7      0.0093121  0.0034175   2.725 0.007027 ** 
## dinfcol.l8   0.1266067  0.0757263   1.672 0.096173 .  
## ddep.l8     -0.0106926  0.0034620  -3.089 0.002310 ** 
## dinfcol.l9  -0.1058791  0.0758295  -1.396 0.164243    
## ddep.l9      0.0065755  0.0035410   1.857 0.064847 .  
## dinfcol.l10  0.0604906  0.0756692   0.799 0.425041    
## ddep.l10    -0.0025046  0.0035172  -0.712 0.477273    
## dinfcol.l11  0.0137712  0.0761681   0.181 0.856715    
## ddep.l11     0.0081791  0.0034852   2.347 0.019954 *  
## dinfcol.l12 -0.3301225  0.0762699  -4.328 2.42e-05 ***
## ddep.l12    -0.0043815  0.0035291  -1.242 0.215927    
## dinfcol.l13  0.1663388  0.0700380   2.375 0.018535 *  
## ddep.l13     0.0025559  0.0038297   0.667 0.505331    
## const        0.0003006  0.0170698   0.018 0.985967    
## dinfeeuu     0.1537803  0.0396354   3.880 0.000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.2509 on 192 degrees of freedom
## Multiple R-Squared: 0.5564,  Adjusted R-squared: 0.4941 
## F-statistic:  8.92 on 27 and 192 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation ddep: 
## ===================================== 
## ddep = dinfcol.l1 + ddep.l1 + dinfcol.l2 + ddep.l2 + dinfcol.l3 + ddep.l3 + dinfcol.l4 + ddep.l4 + dinfcol.l5 + ddep.l5 + dinfcol.l6 + ddep.l6 + dinfcol.l7 + ddep.l7 + dinfcol.l8 + ddep.l8 + dinfcol.l9 + ddep.l9 + dinfcol.l10 + ddep.l10 + dinfcol.l11 + ddep.l11 + dinfcol.l12 + ddep.l12 + dinfcol.l13 + ddep.l13 + const + dinfeeuu 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## dinfcol.l1  -0.6317057  1.2862982  -0.491  0.62391    
## ddep.l1      0.0190786  0.0726720   0.263  0.79319    
## dinfcol.l2  -0.1263818  1.4345492  -0.088  0.92989    
## ddep.l2     -0.0602711  0.0649275  -0.928  0.35443    
## dinfcol.l3  -0.7792393  1.4055971  -0.554  0.57996    
## ddep.l3      0.0068875  0.0645522   0.107  0.91514    
## dinfcol.l4   2.0955911  1.4043602   1.492  0.13729    
## ddep.l4      0.0009237  0.0642012   0.014  0.98854    
## dinfcol.l5  -0.0379583  1.4117828  -0.027  0.97858    
## ddep.l5      0.0475790  0.0646624   0.736  0.46275    
## dinfcol.l6   0.1181141  1.4052348   0.084  0.93310    
## ddep.l6      0.0237484  0.0637802   0.372  0.71004    
## dinfcol.l7  -0.1723914  1.4027883  -0.123  0.90232    
## ddep.l7     -0.0597443  0.0641294  -0.932  0.35270    
## dinfcol.l8   2.4568168  1.4210167   1.729  0.08543 .  
## ddep.l8     -0.0638356  0.0649656  -0.983  0.32704    
## dinfcol.l9   1.1050323  1.4229527   0.777  0.43836    
## ddep.l9      0.0152132  0.0664470   0.229  0.81915    
## dinfcol.l10 -0.3975975  1.4199446  -0.280  0.77977    
## ddep.l10    -0.0028609  0.0660007  -0.043  0.96547    
## dinfcol.l11 -2.3733258  1.4293078  -1.660  0.09845 .  
## ddep.l11     0.0617132  0.0654002   0.944  0.34655    
## dinfcol.l12  0.8665130  1.4312176   0.605  0.54560    
## ddep.l12    -0.4214692  0.0662241  -6.364 1.41e-09 ***
## dinfcol.l13 -2.6467823  1.3142746  -2.014  0.04542 *  
## ddep.l13     0.0468258  0.0718657   0.652  0.51546    
## const        0.0578397  0.3203171   0.181  0.85690    
## dinfeeuu    -1.9600375  0.7437649  -2.635  0.00909 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 4.709 on 192 degrees of freedom
## Multiple R-Squared: 0.3396,  Adjusted R-squared: 0.2468 
## F-statistic: 3.657 on 27 and 192 DF,  p-value: 7.398e-08 
## 
## 
## 
## Covariance matrix of residuals:
##         dinfcol     ddep
## dinfcol 0.06296  0.02799
## ddep    0.02799 22.17028
## 
## Correlation matrix of residuals:
##         dinfcol    ddep
## dinfcol 1.00000 0.02369
## ddep    0.02369 1.00000

Estabilidad del VAR

barplot(roots(var1.est), ylim = c(0,1),
        main = "Estabilidad: Raíces de Polinomio Característico",
        names.arg = as.character(1:2))
abline(h = 1, col = "red", lty = 2, lwd = 2)

barplot(roots(var13.est), ylim = c(0,1),
        main = "Estabilidad: Raíces de Polinomio Característico",
        names.arg = as.character(1:26))
abline(h = 1, col = "red", lty = 2, lwd = 2)

En los dos modelos VAR las raíces son menores que 1, lo que indica estabilidad

Chequeando los residuales de los VAR

La evaluación de los residuales permitirá determinar cual de los modelos se ajusta mejor

plot(var1.est, names = "dinfcol")

plot(var1.est, names = "ddep")

plot(var13.est, names = "dinfcol")

plot(var13.est, names = "ddep")

Al parecer el VAR(13) presenta un mejor ajustes. Miremos pruebas formales de autocorrelación, heteroscedasticidad y normalidad

Autocorrelación

var1.serial1 <- serial.test(var1.est, lags.pt = 20, type = "PT.asymptotic")
var1.serial1
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1.est
## Chi-squared = 162.41, df = 76, p-value = 3.202e-08
var13.serial <- serial.test(var13.est, lags.pt = 20, type = "PT.asymptotic")
var13.serial
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var13.est
## Chi-squared = 52.846, df = 28, p-value = 0.00306

En los dos modelos se rechaza la \(H_0\): no autocorrelación, lo cual indica problemas de autocorrelación

Heteroscedasticidad

var1.hetero <- arch.test(var1.est, lags.multi = 5)
var1.hetero
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object var1.est
## Chi-squared = 93.377, df = 45, p-value = 3.087e-05
var13.hetero <- arch.test(var13.est, lags.multi = 5)
var13.hetero
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object var13.est
## Chi-squared = 62.139, df = 45, p-value = 0.04589

Se observa que se rechaza \(H_0\): homoscedasticidad, lo cual indica problemas de heteroscedasticidad

Normalidad

var1.norm <- normality.test(var1.est)
var1.norm
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var1.est
## Chi-squared = 63.779, df = 4, p-value = 4.651e-13
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var1.est
## Chi-squared = 0.96046, df = 2, p-value = 0.6186
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var1.est
## Chi-squared = 62.819, df = 2, p-value = 2.287e-14
var13.norm <- normality.test(var13.est)
var13.norm
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var13.est
## Chi-squared = 41.83, df = 4, p-value = 1.809e-08
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var13.est
## Chi-squared = 8.4586, df = 2, p-value = 0.01456
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var13.est
## Chi-squared = 33.371, df = 2, p-value = 5.67e-08

Se observa que los residuales de los modelos no son normales (se rechaza \(H_0\): normalidad)

Test de causalidad de Granger

var1.cause.infcol <- causality(var1.est, cause = "ddep")
var1.cause.infcol
## $Granger
## 
##  Granger causality H0: ddep do not Granger-cause dinfcol
## 
## data:  VAR object var1.est
## F-Test = 22.379, df1 = 1, df2 = 456, p-value = 2.99e-06
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: ddep and dinfcol
## 
## data:  VAR object var1.est
## Chi-squared = 1.0761, df = 1, p-value = 0.2996
var13.cause.infcol <- causality(var13.est, cause = "ddep")
var13.cause.infcol
## $Granger
## 
##  Granger causality H0: ddep do not Granger-cause dinfcol
## 
## data:  VAR object var13.est
## F-Test = 2.8431, df1 = 13, df2 = 384, p-value = 0.0006401
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: ddep and dinfcol
## 
## data:  VAR object var13.est
## Chi-squared = 0.1234, df = 1, p-value = 0.7254

El test de Granger muestra que la depreciación causa Granger la inflación

Predicción

Para la predicción de la inflación para Colombia y la depreciación, se deben asignar valores futuros para la variable exogena inflación de los Estados Unidos y luego diferenciarlos. Se toman los datos de inflación en los Estados Unidos de junio de 2022 a enero de 2023 (se toman 8 meses que al diferenciar nos da 7 meses), lo cual daría para hacer 7 meses de predicción.

infeeuu_f <- as.matrix(data.frame(dinfeeuu = c(8.93298689, 8.413182026, 8.227361014, 8.214853957, 7.762492677,
                                               7.135348085, 6.444940492, 6.347156218)))
dinfeeuu_f <- diff(infeeuu_f,1)
var1.predic <- predict(var1.est, n.ahead = 7, ci = 0.95, dumvar = dinfeeuu_f)
plot(var1.predic, names = "dinfcol")

plot(var1.predic, names = "ddep")

data.predict1 <- data.frame(date = as.yearmon(seq(as.Date("2003-01-01"), as.Date("2023-01-01"), by = "month"), "%m-%Y"),
                            predicho = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var1.predic$fcst$dinfcol[1:7,1]),
                            ic_l = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var1.predic$fcst$dinfcol[1:7,2]),
                            ic_u = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var1.predic$fcst$dinfcol[1:7,3]),
                            actual = c(data$dinfcol, rep(NA, 7)))

ggplot(data.predict1[data.predict1$date>"ene. 2020",]) + 
  geom_line(aes(x = date, y = actual, color = "Observada"), size = 0.8) + 
  geom_line(aes(x = date, y = predicho, color = "Predicha"), size = 0.8) + 
  geom_ribbon(aes(x = date, y = predicho, ymin = ic_l, ymax = ic_u, fill="IC"), alpha = 0.1) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(title = "Predicción del D(inflación) de Colombia VAR(1)", x = "Mes", y ="", fill = "") + 
  scale_x_continuous(n.breaks = 8)

var13.predic <- predict(var13.est, n.ahead = 7, ci = 0.95, dumvar = dinfeeuu_f)
plot(var13.predic, names = "dinfcol")

plot(var13.predic, names = "ddep")

data.predict13 <- data.frame(date = as.yearmon(seq(as.Date("2003-01-01"), as.Date("2023-01-01"), by = "month"), "%m-%Y"),
                            predicho = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var13.predic$fcst$dinfcol[1:7,1]),
                            ic_l = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var13.predic$fcst$dinfcol[1:7,2]),
                            ic_u = c(rep(NA, length(data$dinfcol) - 1), data$dinfcol[length(data$dinfcol)],
                                               var13.predic$fcst$dinfcol[1:7,3]),
                            actual = c(data$dinfcol, rep(NA, 7)))

ggplot(data.predict13[data.predict13$date>"ene. 2020",]) + 
  geom_line(aes(x = date, y = actual, color = "Observada"), size = 0.8) + 
  geom_line(aes(x = date, y = predicho, color = "Predicha"), size = 0.8) + 
  geom_ribbon(aes(x = date, y = predicho, ymin = ic_l, ymax = ic_u, fill="IC"), alpha = 0.1) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(title = "Predicción del D(inflación) de Colombia VAR(13)", x = "Mes", y ="", fill = "") + 
  scale_x_continuous(n.breaks = 8)

Lo anterior fue la predicción de la serie de inflación diferenciada. Ahora, para hacer la predicción de la serie original se procede como sigue:

\[\Delta x_t = x_t - x_{t-1} \] \[\Delta x_{t+1} = x_{t+1} - x_{t} \]

Así que, para deducir el valor predicho en niveles seria: \[x_{t+1} = \Delta x_{t+1} + x_{t}\]

Con los datos predichos de los modelos VAR en diferencias, se puede hacer el cálculo de la prediccion de la variable original de inflación

predict1_infcol <- data.frame(infcol_p = NA)

for (i in 1:7) {
predict1_infcol <- data.frame(infcol_p = c(
  data.predict1$predicho[235]+data$infcol[234],
  data.predict1$predicho[236]+predict1_infcol$infcol_p[1],
  data.predict1$predicho[237]+predict1_infcol$infcol_p[2],
  data.predict1$predicho[238]+predict1_infcol$infcol_p[3],
  data.predict1$predicho[239]+predict1_infcol$infcol_p[4],
  data.predict1$predicho[240]+predict1_infcol$infcol_p[5],
  data.predict1$predicho[241]+predict1_infcol$infcol_p[6]))
}

data.predict1infcol <- data.frame(date = as.yearmon(seq(as.Date("2003-01-01"), as.Date("2023-01-01"), by = "month"), "%m-%Y"),
                                  predicho = c(rep(NA, length(data$infcol) - 1), data$infcol[length(data$infcol)], predict1_infcol$infcol_p),
                                  actual = c(data$infcol, rep(NA, 7)))

ggplot(data.predict1infcol) + 
  geom_line(aes(x = date, y = actual, color = "Observada"), size = 0.8) + 
  geom_line(aes(x = date, y = predicho, color = "Predicha"), size = 0.8) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  labs(title = "Predicción de la inflación de Colombia VAR(1)", x = "Mes", y ="", fill = "") + 
  scale_x_continuous(n.breaks = 8) + scale_y_continuous(n.breaks = 6)

predict13_infcol <- data.frame(infcol_p = NA)

for (i in 1:7) {
predict13_infcol <- data.frame(infcol_p = c(
  data.predict13$predicho[235]+data$infcol[234],
  data.predict13$predicho[236]+predict13_infcol$infcol_p[1],
  data.predict13$predicho[237]+predict13_infcol$infcol_p[2],
  data.predict13$predicho[238]+predict13_infcol$infcol_p[3],
  data.predict13$predicho[239]+predict13_infcol$infcol_p[4],
  data.predict13$predicho[240]+predict13_infcol$infcol_p[5],
  data.predict13$predicho[241]+predict13_infcol$infcol_p[6]))
}

data.predict13infcol <- data.frame(date = as.yearmon(seq(as.Date("2003-01-01"), as.Date("2023-01-01"), by = "month"), "%m-%Y"),
                                  predicho = c(rep(NA, length(data$infcol) - 1), data$infcol[length(data$infcol)], predict13_infcol$infcol_p),
                                  actual = c(data$infcol, rep(NA, 7)))

ggplot(data.predict13infcol) + 
  geom_line(aes(x = date, y = actual, color = "Observada"), size = 0.8) + 
  geom_line(aes(x = date, y = predicho, color = "Predicha"), size = 0.8) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  labs(title = "Predicción de la inflación de Colombia VAR(13)", x = "Mes", y ="", fill = "") + 
  scale_x_continuous(n.breaks = 8) + scale_y_continuous(n.breaks = 6)

Funciones impulso respuesta

irf1.dinfcol <- irf(var1.est, impulse = "ddep", response = "dinfcol", 
                  n.ahead = 10, ortho = F, boot = TRUE, runs = 1000,
                  seed = 0123456789)

data.irf1 <- data.frame(time = seq(1,11,1),
                       irf = irf1.dinfcol$irf$ddep[1:11,1],
                       lower = irf1.dinfcol$Lower$ddep[1:11,1],
                       upper = irf1.dinfcol$Upper$ddep[1:11,1])

ggplot(data.irf1) + 
  geom_line(aes(x = time, y = irf, color = "FIR"), size = 0.8) + 
  geom_ribbon(aes(x = time, y = irf, ymin = lower, ymax = upper, fill="IC"), alpha = 0.1) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("FIR" = "darkblue")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(title = "Shock de la depreciación sobre la inflación VAR(1)", x = "Mes", y ="Inflación de Colombia", fill = "", caption ="Intervalos de confianza con 1000 bootstrap") + 
  scale_x_continuous(n.breaks = 11) +
  geom_hline(yintercept=0, color = "red")

irf13.dinfcol <- irf(var13.est, impulse = "ddep", response = "dinfcol", 
                  n.ahead = 10, ortho = F, boot = TRUE, runs = 1000,
                  seed = 0123456789)

data.irf13 <- data.frame(time = seq(1,11,1),
                       irf = irf13.dinfcol$irf$ddep[1:11,1],
                       lower = irf13.dinfcol$Lower$ddep[1:11,1],
                       upper = irf13.dinfcol$Upper$ddep[1:11,1])

ggplot(data.irf13) + 
  geom_line(aes(x = time, y = irf, color = "FIR"), size = 0.8) + 
  geom_ribbon(aes(x = time, y = irf, ymin = lower, ymax = upper, fill="IC"), alpha = 0.1) +
  theme(legend.text = element_text(size = 10), text = element_text(size=12), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("FIR" = "darkblue")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(title = "Shock de la depreciación sobre la inflación VAR(13)", x = "Mes", y ="Inflación de Colombia", fill = "", caption ="Intervalos de confianza con 1000 bootstrap") + 
  scale_x_continuous(n.breaks = 11) +
  geom_hline(yintercept=0, color = "red")

Referencias

Posada, CE y Londoño-Sierra, L. (2022). “Retos para la autoridad monetaria en el próximo cuatrenio”. Grupo de Coyuntura Económica. Universidad EAFIT, Medellín-Colombia. https://coyunturaeconomicaeafit.wordpress.com/2022/08/25/retos-para-la-autoridad-monetaria-en-el-proximo-cuatrienio/